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ABSTRACT 


The  objective  is  to  obtain  and  compare  interpolation 
methods  and  contour  plot  techniques ,  based  on  data  from 
sites  whose  locations  are  irregular  or  scattered.  Data 
from  the  1980  Fugro  report  ^Ref.  lj  on  MX  valley  soil 
samples  serves  as  a  test  bed. 

Two  groups  of  data,  seismic  p-wave  velocities  and  sur¬ 
face  soil  depth  were  studied  using  six  different  interpola 
tion  methods  and  three  different  contour  plot  techniques. 
Comparison  of  the  interpolation  methods  and  of  the  contour 
plot  techniques  are  made. 
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INTRODUCTION 


I . 

Detailed  investigations  of  seismic  p-wave  velocities  and 
depth  of  the  surface  soil  layer  were  performed  in  the  Ralston 
Valley,  Nevada,  by  Fugro ,  Inc. [Ref.  ij  ,  for  the  U.S.  Army 
Engineer  Waterways  Experiment  Station  (W.E.S.)  as  part  of  the 
preliminary  work  on  the  MX  missile  system  deployment.  Six¬ 
teen  locations  were  chosen  in  the  Ralston  Valley  by  W.E.S. 
for  their  study  organized  into  four  sets  of  sites  (coded  RA, 
RB , RC , RD)  located  in  each  of  the  four  predominant  surficial 
soil  types  (coded  5I,5Y,U,4U),  and  an  additional  site  marked 
RU2  identified  in  the  Valley. 

The  data  source  for  identified  sites  Qtef.  2J  gives  the 
(X, Y)  coordinates  (one  mile  to  the  inch)  of  the  boring 
locations  relative  to  the  map's  lower  left  grid  mark,  the 
seismic  p-wave  velocities  (MPS  -  meters  per  second)  of  the 
surface  layer,  and  the  depth  of  that  layer  (meters)  for  the 
2B  line;  see  Appendix  A.  The  site  identifiers  and  locations 
for  the  Ralston  Valley  sampling  plan  can  be  found  in  Appendix 
B.  The  immediate  concern  is  to  develop  contours  describing 
the  depth  of  the  top  layer,  and  of  the  compression  wave  speed 
in  the  top  layer,  in  the  entire  Valley. 

Based  on  these  limited  data,  the  mathematical  problem  is 
that  of  constructing  a  smooth  bivariate  function  F(X,Y)  which 

Y  )  =  Fk  for  K  =  1 , . . . ,N. 


takes  on  specified  values  FfX^f 
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Franke  ([Ref.  3,4,5,6j  has  made  a  rather  extensive  general 
study  of  the  modern  interpolation  methods.  Most  of  the  methods 
require  the  user  to  specify  a  parameter  (or  several) ,  but 
these  parameters  can  be  determined  by  computational  experi¬ 
ment  ^Ref .  4 J  .  The  interpolation  methods  selected  for 
specific  inclusion  in  the  current  work  are  a  subset  of  these 
methods,  see  Table  1.  Furthermore,  contour  plotting,  specific¬ 
ally  non-IMSL  subroutine  CONTUR,  subroutine  CONISD  and  sub¬ 
routine  PLT3D1,  were  exercised  for  comparison.  Basic  FORTRAN 
programs  for  the  first  five  interpolation  methods  were  obtained 
from  Professor  R.  Franke  [[Ref.  3, 4,5,6  J. 

The  various  interpolation  methods  were  used  to  supply 
function  values  at  a  gridwork  of  reference  points.  In  order 
to  compare  the  methods,  the  mean  and  standard  deviation  of 
the  values  were  produced  at  each  grid  point.  Also  it  is 
important  to  identify  which  methods  generates  sharply  different 
interpolation  values  compared  to  the  other  methods.  For  this 
purpose,  the  same  FORTRAN  program  generates  the  absolute  dif¬ 
ferences  from  the  mean  values  for  each  method. 

FORTRAN  programs  for  all  of  the  methods  can  be  obtained 
by  writing  the  advisor. 

In  the  next  chapter  we  will  describe  a  number  of  inter¬ 
polation  methods.  Chapter  3  describes  employment  of  the 
Ralston  Valley  data  as  a  way  of  comparing  these  methods. 

Chapter  4  discusses  contour  plot  techniques ,  and  the  con¬ 
clusions  from  using  these  procedures  are  given  in  Chapter  5. 
Graphical  results  are  presented  in  a  series  of  appendices. 


TABLE  1 


INTERPOLATION  METHODS 

1)  Inverse  Distance  Weighted  Methods 

A)  Shepard's  Method 

B)  Modified  Shepard's  Method 

C)  Modified  Linear  Shepard's  Method 

D)  Modified  Shepard's  Method  Boolean  Sum  Plane 

E)  Modified  McLain  Method 

F)  Quadratic  Shepard's  Method 

*G)  Modified  Quadratic  Shepard's  Method 
H)  McLain  Method 

2)  Franke's  Methods  (Rectangle  based  blending  methods) 

A)  Franke's  Method  (Mode  One) 

B)  Franke's  Method  (Mode  Three) 

*C)  Franke's  Method  (Thin  plate  local  functions) 

3)  Triangle  Based  Blending  Methods 

A)  Nielson-Franke  Linear  Triangle  Method 
*B)  Nielson-Franke  Quadratic  Triangle  Method 

4)  Finite  Element  Based  Methods 

A)  Akima's  Method 

B)  Akima's  Method  -  Modification  Two 

C)  Akima's  Method  -  Modification  One 

D)  Akima's  Method  -  Modification  Three 
*E)  Nielson's  Minimum  Norm  Network 

F)  Lawson's  Method 

5)  Foley's  Method 

A)  Generalized  Newton  Interpolant 

B)  TF  Delta  Sum  Berstein  Interpolant 

C)  Iterated  Delta  Sums :  TF  Delta  Sum  Bicubic  Spline 
*D)  Iterated  Delta  Sums:  A  Shepard's  Method  Delta  Sum 

Bic.  Sp. 

6)  Nodal  Basis  Function  Type  Methods 
A)  Rotated  Gaussian 

*B)  Hardy's  Multiquadric 

C)  Hardy's  Reciprocal  Multiquadric 

D)  Duchon's  Radial  Cubic  Method 

E)  Duchon's  Thin  Plate  Function 

F)  Rotated  B-Spline 


*  Methods  chosen  for  comparison  study. 


II.  INTERPOLATION  METHODS 


The  problem  being  addressed  here  is  that  of  constructing 
a  smooth  bivariate  function,  F(X,Y) ,  which  takes  on  specified 
values  F(Xk,Y^)  =  F^,  k=l,...,N  based  on  point  data  from 
scattered,  irregular  locations.  The  bivariate  function, 

F (X ,Y)  is  smooth,  {i.e.  ,  it  has  at  least  continuous  first 
partial  derivatives),  and  the  points  (X^Y^)  are  * scattered ’ > 
that  is,  not  necessarily  conforming  to  some  regular  pattern. 

Some  interpolation  methods,  called  'global',  are  sensitive 
to  all  the  data  points :  the  addition  or  deletion  of  a  data 
point,  or  a  change  of  one  of  the  coordinates  of  a  data  point, 
will  propogate  throughout  the  entire  domain  of  definition. 

Other  interpolation  methods  are  called  'local',  that  is,  insen¬ 
sitive  to  changes  in  the  data  points  which  are  sufficiently 
distant  from  the  location  at  hand. 

Franke,  R.  £Ref.  3,4,5,6j  has  made  an  extensive  general 
study  of  new  interpolation  methods.  He  produced  a  classifi¬ 
cation  of  the  interpolation  methods  [Ref.  4],  as  shown  in 
Table  1.  Under  his  six  main  types,  there  are  a  total  of  29 
different  methods.  All  the  methods  appearing  in  the  same 
group  are  either  local  or  global  and  use  the  same  basic  idea; 
several  are  modifications  of  the  previously  described  ones. 

Six  methods ,  one  from  each  of  the  six  main  types ,  are 


described  below. 


A.  MODIFIED  QUADRATIC  SHEPARD'S  METHOD 


The  modified  quadratic  Shepard's  method  was  derived 
from  the  inverse  distance  weighted  method.  All  methods  of 
this  type  which  are  considered  may  be  viewed  as  generaliza¬ 
tions  of  Shepard's  method.  The  basic  Shepard's  method  is 

N  N 

F(X,Y)  =  Z  W.  (X  ,Y)  .  Q.  (X,Y)  /  Z  W.(X,Y) 
k=l  *  k=l 

“U 

where  w,(X,Y)  =  d,  .  Typically  p  =  2 ,  although  other  values 
k  k 

may  be  used  Here  d2  =  ( (X-K^) 2  +  (Y-Y^) 2) . 

In  the  modified  quadratic  Shepard's  method  the  weights 
for  obtaining  the  nodal  functions  (quadratics)  are  taken  as 

2 

and  R  =  1/2  l/N/N.'D 

q  q 

where  D  is  the  diameter  of  the  point  set,  and  a  nominal  value 
for  Ng,  determined  by  computational  experiment  £ Ref .  4],  is 
N  =  18.  The  diameter  of  the  point  set  is  defined 

q 

2  f  2  2 

D  =  MAX  (X.  -  X.)  +  (Y.  -  Y  ) 

I,J  L  1  3  1  ^  J 

The  nodal  functions  Q  (X,Y)  are  defined  as 

Jv 


Wfc(X,Y)  « 


(Rq  “  V 


(R, 


V 
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rWi'Ai«ii! 


2 


Q.(X,Y)  =  F  +  a  (X-X  )  +  a,  (Y-Y,  }  +  a  (X-X  ) 
k  X  k2  k  k3  X  k4  k 


+  a  (X-X  ) •  (Y-Y  )  +  a  (Y-Y  ) 
k5  k  k  X6  k 


and  for  the  coefficients  one  must  solve  the  least  squares 


problem. 


JL  W*  [Flc  *  WVV  +  *MlW 


a,  .  ,  j— 2  , .  • .  ,  6  i£k 

+ 


\4(W  +akS(W-(W 


+  ^‘VV  -  Fi] 


Complete  details  are  given  in  £  Ref.  4,5  J. 


B.  LOCAL  THIN  PLATE  SPLINES  METHOD 


This  is  known  as  Franke's  method  £  Ref.  4  J.  The  idea  is 


to  represent  the  interpolation  function  as 


F  (X ,Y)  =  Z  W.  (X ,Y)  -Q  .  (X  ,Y) 
i ,  i  V 


where  the  W.  .  (X,Y)  are  local  weight  functions  I"  Ref.  6l.  Here 

we  choose  the  W  so  that  W. .  =  1,  and  the  local  approximations 
i,j  **3 

Q, j (X,Y)  interpolate  at  points  where  the  corresponding  weight 
-V 

function  is  non-2ero.  The  weight  functions  are  products  of 
piecewise  hermite  cubicsJTRef.  fij. 
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The  local  approximations  for  this  method  are  taken  to  be 


the  thin  plate  splines  and  have  the  form 

2 

Q.- (X,Y)  =  Z  -LOG  d.  +  a  +  b-X  +  C-Y 

^  k£l  ^  k  k 

2  2  2 
where  d£  =  ( (X-Xk)  +  (Y-Yk)  )  . 

And  I  is  the  set  of  indices  k  for  which  (X^ ,Yk ,F^)  is  a  point 
to  be  interpolated  by  Qj_j(X,Y) ,  and  is  defined  by 

is  to  take  the  value  at  ( 

The  coefficients  A^,  and  a,  b  and  c  are  determined  by  a  certain 
linear  system  of  equations,  (see  [[Ref.  6j)  . 

C.  NIELSON-FRANKE  QUADRATIC  TRIANGLE  METHOD 

All  methods  listed  under  triangle  based  blending  methods 
J^Ref.  4j  are  conceptually  similar  to  the  inverse  distance 
weighted  methods.  The  interpolation  function  for  the  quad¬ 
ratic  triangle  method  is 

F  (X,Y)  -  W^X/Y) -Q  (X,Y)  +  Wj  (X,Y)  *Qj  (X,Y)  +  Wk (X,Y) -Qk (X ,Y) 

The  first  step  in  the  present  method  is  to  partition  the  plane 
into  triangles  by  connecting  neighboring  data  points  based 
upon  the  min-max  angle  criterion  as  described  by  Franke  in 
£*Ref.  5j.  The  quadratic  triangle  method  uses  the  inverse 
distance  weighted  quadratic  Q^(X,Y)  which  is  also  used  in  the 
modified  quadratic  Shepard's  method,  see  section  2. A. 


A  significant  difference  from  the  inverse  distance  weighted 
method  is  that  the  weight  functions  for  the  triangle  based 
blending  method  are  based  on  a  triangulation  of  the  convex 

'V 

tive  to  the  triangulation  technique.  If  a  triangulation  is  in 
existence  for  other  purposes ,  it  can  be  used. 

In  short,  the  method  is  to  define  the  nodal  functions 
,k=*l , . . .  ,N ,  as  in  the  modified  quadratic  Shepard's  method, 
form  a  triangulation  of  the  points  ,  determine  the  weight 
functions  [^Ref.  5  J  and  compute  interpolant  F(X,Y). 

D.  NIELSON'S  MINIMUM  NORM  NETWORK  METHOD 

Nielson's  minimum  norm  network  method  ^Ref.  4,5j,  can 
be  found  under  the  finite  element  based  methods  in  Table  1. 

Like  the  Nielson-Franke  quadratic  triangle  method,  this  global 
method  is  based  upon  a  tri angulation. 

The  method  consists  of  three  separate  steps;  triangulation 
curve  network,  and  blending  [^Ref.  7J.  The  points 
(X  ,Y.)  ,k=l , . . .  ,N  are  used  as  the  vertices  of  a  triangulation 

K  & 

of  the  convex  hull  of  these  points.  Alternatively,  if  a  tri- 
aggulation  already  exists,  it  can  be  used.  The  curve  network 
step  involves  the  solution  of  a  certain  minimum  norm  problem 
and  requires  first  partial  derivatives  in  its  discretized 
form.  These  are  obtained  by  assuming  a  cubic  variation  along 
each  edge  in  the  tri angulation  and  minimizing  the  integral  of 
the  second  derivative  squared  ^Ref.  7J. 


hull  of  the  point  set  ^ (X^ 


} 


The  method  is  not  sensi- 


19 


This  method  does  not  provide  extrapolation  outside  the  con¬ 
vex  hull,  but  the  triangular  blending  method  can  be  used  to 
extend  the  curve  network  to  the  entire  hull.  In  a  triangular 
exterior  region,  the  function  is  taken  to  be  the  linear 
function  determined  by  the  value  and  slopes  at  the  vertex 
[Ref.  7]. 


E.  GENERALIZED  NEWTON  DELTA  SUM  BICUBIC  SPLINE  METHOD 

Foley's  methods  [Ref.  4J  involve  several  ideas.  The 
interpolant  is  taken  as  either  the  generalized  Newton  inter- 
polant,  or  a  form  of  Shepard's  method.  The  use  of  generalized 
Newton  type  interpolant  is  involved  in  them  prominently, 
because  the  best  performance  for  various  test  functions  is 
provided  by  the  interated  delta  sum  methods  using  the  generalized 
Newton  interpolant  with  natural  bicubic  splines. 

The  generalized  Newton  interpolant  takes  the  form 

N 

T  (X ,Y)  =  l  a,  -Wk(X,Y)  f 
N  k=l  k 


where 


ak  = 


f  -  T  (x  ,y  ) 
k  k-1  k  k 


vw 


and  W  ( X , Y )  has  the  property  W  (X. ,Y. )  =  0  i=l , 2 , . . . ,k-l  . 
k  k  1  1 

This  function  is  dependent  upon  the  order  of  the  points,  and 
so  Foley's  scheme  involves  an  ordering  process  [Ref.  4J. 
After  constructing  a  bicubic  spline  interpolant  for  the  grid 
points,  one  adds  a  generalized  Newton  interpolant  for  the 
difference  between  data  and  the  bicubic  splines,  obtaining 
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an  interpolant.  This  process  is  termed  a  'delta  sum'  by 
Foley  ^Ref.  4J. 


F.  HARDY'S  MULTIQUADRIC  METHOD 


Hardy's  multiquadric  method  [Ref.  8,9^  can  be  thought 
of  as  being  under  the  global  basis  function  type.  A  multi¬ 
quadric  surface  can  be  represented  by 


N  2  2 

Z  c.  (X.-X. )  +  (Y  -Y  )  +  c 

j=l  1  3  1  3  i 


—  Z^ /  i— 1  / . •  .  ,N 


The  given  problem  data  provide  a  set  of  cartesian  coordinates 
on  the  surface  ranging  from  XjyYj^Z^  to  X^Y^z^-  the  quad¬ 
ratic  term  coefficients  to  are  unknown.  The  coordinates 
of  the  N  data  points  can  be  substituted  into  the  definition 
of  the  surface  and  only  the  coefficients  C  are  left  as 

j 

unknown.  A  problem  here  is  that  of  solving  a  system  of  N 
linear  equations  with  N  unknowns .  A  nominal  value  for  C  was 
determined  by  computational  experiment  and  the  value  C  =  0 
was  chosen. 


If  a  vector  of  unknowns  is  to  be 
used,  define 


then  each  known  element  becomes 


t 


I 
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III.  COMPARISON  OF  INTERPOLATION  METHODS  USING 
RALSTON  VALLEY  DATA 

Contour  construction  for  seismic  p-wave  velocities  and 
depth  data  can  be  seen  in  Appendix  C  and  in  Appendix  F , 
respectively,  for  the  various  interpolation  methods  and  with 
the  three  different  contour  plot  techniques.  Each  method 
appears  in  the  order  given  in  Table  1.  Each  method  is  dis¬ 
played  using  first  the  subroutine  CONTUR,  then  the  subroutine 
CON ISP,  and  lastly  the  subroutine  PLT3D1 .  Subroutine  CONTUR 
cannot  outline  the  Ralston  Valley,  but  gives  contour  con¬ 
structions  in  rectangular  form  with  the  x-axis  ranging  from 
one  to  seventeen  miles  and  the  y-axis  from  one  to  twenty  five 
miles.  The  subroutine  CONISD  can  outline  the  Ralston  Valley 
with  boundary  lines.  All  three  dimensional  pictures  were 
generated  by  subroutine  PLT3D1. 

Appendix  D  shows  the  mean  *nd  standard  deviation  of  the 
various  contour  constructions  for  the  seismic  p-wave  velo¬ 
cities.  Numerical  values  for  the  mean  and  standard  deviation 
at  each  grid  point  come  from  six  different  interpolation 
methods.  Appendix  G  is  similar  to  Appendix  D,  but  it  is  for 
the  depth  of  the  surface  soil  layer.  Convex  contour  lines 
are  formed  for  the  standard  deviations  values  as  can  be  seen 
in  Appendices  D  and  G.  Small  standard  deviation  values 
appear  in  data  point  region,  and  become  larger  as  one  moves 
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toward  the  edges.  In  other  words,  more  accurate  prediction 
occurs  in  data  point  region  and  less  accurate  prediction  out¬ 
side  this  region. 

Deviations  from  mean  values  for  all  interpolation  methods 
can  be  seen  in  Appendix  E  for  seismic  p-wave  velocities  and 
in  Appendix  H  for  depth  values.  The  1,  10,  50  and  100  contours 
show  the  absolute  differences  between  the  given  method  and 
mean  values  for  seismic  p-wave  velocities.  The  0.1,  0.5,  1.0 
contours  for  depth  data  display  quite  convex  regions  for  all 
methods.  Clearly,  it  can  be  seen  that; there  are  no  signifi¬ 
cant  differences  between  methods  in  terms  of  deviations  from 
mean  values  for  grid  points  especially  in  data  region.  This 
result  shows  that  in  the  data  points  region  all  methods  work 
with  almost  same  accuracy  and  outside  the  region  there  are 
too  great  differences.  There  is  no  way  to  compare  them  out¬ 
side  the  region. 


IV.  CONTOUR  PLOT  TECHNIQUES 


The  interpolation  values  provided  by  any  of  the  methods 
can  be  seen  in  the  computer  outputs ,  but  their  interpretation 
is  highly  difficult  without  contour  construction.  For  inter¬ 
pretation  purposes,  computer  programs  for  drawing  plots  on 
an  off-line  plotter  are  available  at  the  Naval  Postgraduate 
School  computer  center. 

The  subprograms  CONTUR,  CONISD  and  PLT3D1  serve  the  com¬ 
mon  purpose  but  have  different  capabilities. 

All  the  subprograms  have  the  same  limitation  in  that  they 
use  a  linear  interpolation  process  to  find  the  values  of  the 
points  along  contour  segments.  Also,  there  is  a  limit  in  the 
resolution  of  the  off-line  plotter,  namely,  .005  inch  in  both 
the  X  and  Y  directions  for  the  versatec  plotter. 

A.  SUBROUTINE  CONTUR 

Given  a  matrix  of  numerical  values  of  Z  =  (X,Y) ,  the 
program  CONTUR  generates  a  contour  graph  on  which  one  or 
more  contour  levels  are  drawn.  This  non-IMSL  subroutine 
uses  a  two-dimensional  array  and  turns  out  a  two-dimensional 
off-line  plot.  Interior  and  exterior  contour  segments  and 
labels  (if  requested)  for  exterior  contour  segments  can  be 
seen.  Labeling  of  the  interior  segments  represent  local 
maxima. 
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Advantages : 

1)  The  subroutine  has  a  minimum  number  of  arguments  (10) 
compared  with  other  plot  techniques. 

2)  A  rather  big  two-dimensional  array  can  be  used,  up  to 
1797  data  points. 

3)  It  takes  approximately  one  fiftieth  of  the  CPU  time  that 
would  be  required  by  subroutine  CONISD. 

Disadvantages  ; 

1)  Irregularly  spaced  data  cannot  be  contoured. 

2)  Irregular  boundary  capability  is  not  available,  and  there 
is  no  option  for  a  'cut-out'  area. 

3)  The  width  and  the  height  of  the  contour  graph  must  be 
specified  as  integers,  and  are  not  adjustable  scale 
parameters . 

4)  The  scale  values  (if  requested)  one  inch  apart  on  the 
exterior  frame  of  the  contour  graph  increase  from  north 
to  south.  This  is  not  oriented  properly  for  our  use, 
which  is  a  lower  left  grid  oriented  system. 

5)  Limited  labeling.  The  x-axis  and  y-axis  cannot  be 
labeled. 

6)  There  is  no  option  for  smoothing  the  contour  lines. 

B.  SUBROUTINE  CONISD 

This  non-lMSL  subroutine  produces  a  drawn  contour  map  on 

an  off-line  plotter  for  given  irregularly  spaced  data  points. 

Each  data  point  is  a  triad  of  X,Y  and  Z  values  where 

Z  =  (X,Y) .  There  may  be  one  or  more  'cut-out  areas'  which 


are  not  contour. 
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Advantages : 

1)  Irregularly  spaced  data  can  be  drawn. 

2)  Irregular  boundary  capability  is  available  for  one  or 
more  'cut-out  areas'. 

3)  The  scale  for  the  x-axis  and  y-axis  are  defined  as  real 
numbers,  i.e.,  adjustable  scale. 

4)  X-axis  and  y-axis  can  be  labeled. 

5)  There  is  an  option  for  smoothing  the  contour  lines. 
Disadvantages : 

1)  It  takes  much  more  CPU  time  than  the  other  plot  techniques. 
C.  SUBROUTINE  PLT3D1 

A  subroutine  PLT3D1  produces  three-dimensional  perspec¬ 
tive  or  isometric  plots  on  an  off-line  plotter  for  given  a 
two-dimensional  matrix  of  numerical  values  of  Z  -  (X,Y) . 
Labeling  and  scaling  are  not  available  for  this  reason,  it 
can  be  used  only  for  general  purposes. 

A  memorandum  about  this  subprogram  (jtef.  12  J  can  be 
requested  from  W.R.  Church  computer  center.  This  memorandum 
includes  a  listing  of  the  computer  program,  the  algorithm 
that  is  used  by  the  program,  and  the  axes  orientation,  rota¬ 
tion  and  projection  are  clarified  by  illustrations. 

Advantages : 

1)  This  routine  gives  a  nice  overview  to  given  area  by 
Z  *  (X  ,Y)  . 

2)  The  vantage  point  can  be  changed. 

3)  The  scale  is  adjustable. 
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Disadvantages : 


1)  There  is  no  option  for  'cut-out  areas'. 

2)  There  is  no  labeling. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Seismic  p-wave  velocities  or  depth  values  for  any  grid 
point  in  Ralston  Valley  can  be  found  from  computer  outputs  or 
contour  constructions  for  any  particular  interpolation  method 
or  for  the  mean  of  the  six  different  metnods .  The  standard 
deviations  for  each  grid  point  are  also  available.  The  means 
and  standard  deviations  may  be  useful  for  the  comparison  of 
methods . 

The  computer  outputs  for  standard  deviations  in  Appendix 
E  and  H  show  that  accuracy  is  high  in  the  data  point  region 
and  at  some  undecided  distance  from  this  region;  but  this 

! 

distance  cannot  be  measured.  In  the  data  point  region  all 
interpolation  methods  gives  essentially  the  same  result, 
regardless  of  whether  they  are  global  or  local. 

Outside  the  data  point  region  all  methods  generate  rather 
different  values  and  there  is  no  way  to  decide  which  one  is 
close  to  true  values.  If  the  point  of  interest  is  in  the  data 
point  region  any  of  the  interpolation  methods  can  be  used. 

According  to  Appendix  C  and  P;  none  of  the  methods  can 
find  a  maximum  or  minimum  point  outside  the  region.  Table  2 
gives  approximate  CPU  times  for  all  of  the  studied  inter¬ 
polation  methods.  If  minimal  computer  time  is  required  for  a 
given  problem,  the  triangle  based  quadratic  method  is  suitable. 
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APPENDIX  A 


TABLE  3 

DATA  SELECTED  FOR  TEST  BED 


Side 

Designation 

Coordinates 

X  Y 

P-Wave  Vel. 
NPS 

Depth 

M. 

1. 

RA5I 

9.25 

30.25 

561 

3.0 

2. 

RB5I 

5.69 

25.5 

366 

1.0 

3. 

RC5I 

4.75 

20.81 

279 

0.9 

4. 

RD5I 

5.75 

14.75 

347 

1.5 

5. 

RA5Y 

5.59 

26.25 

341 

1.5 

6. 

RB5Y 

11.75 

24.19 

366 

0.8 

7. 

RC5Y 

7.06 

20.0 

323 

1.2 

8. 

RD5Y 

7.31 

10.31 

329 

0.8 

9. 

RAU 

7.31 

28. 19 

363 

1.2 

10. 

RBU 

9.25 

19.31 

335 

1.6 

11. 

RCU 

8.75 

11.25 

317 

1.4 

12. 

RDU 

17.44 

10.94 

378 

1.5 

13. 

RA4U 

10.12 

23.44 

341 

1.5 

14. 

RB4U 

12.25 

18.56 

372 

1.5 

15. 

RC4U 

11.44 

14.38 

372 

1.5 

16. 

RD4U 

15.44 

12.50 

381 

0.5 

17. 

RU2 

6.00 

27.19 

357 

1.2 
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